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This paper addresses the issues of heat transfer in oscillatory flow conditions, which are typically found 
in thermoacoustic devices. The analysis presented concerns processes taking place in the individual 
“channels” of the parallel-plate heat exchangers (HX), and is a mixture of experimental and numerical 
approaches. In the experimental part, the paper describes the design of experimental apparatus to study 
the thermal-fluid processes controlling heat transfer in thermoacoustic heat exchangers on the micro¬ 
scale of the individual channels. Planar Laser Induced Fluorescence (PLIF) and Particle Image Velocim- 
etry (PIV) techniques are applied to obtain spatially and temporally resolved temperature and velocity 
fields within the HX channels. The temperature fields allow obtaining the local and global, phase- 
dependent heat transfer rates and Nusselt numbers, and their dependence on the Reynolds number of 
the oscillating flow. The numerical part of the paper deals with the implementation of CFD modelling 
capabilities to capture the physics of thermal-fluid processes in the micro-scale and to validate the 
models against the experimental data. A two-dimensional low Mach number computational model is 
implemented to analyse the time-averaged temperature field and heat transfer rates in a representative 
domain of the HXs. These are derived by integrating the thermoacoustic equations of the standard linear 
theory into a numerical calculus scheme based on the energy balance. The comparisons between the 
experimental and numerical results in terms of temperature and heat transfer distributions suggest that 
the optimal performance of heat exchangers can be achieved when the gas displacement amplitude is 
close to the length of hot and cold heat exchanger. Heat transfer coefficients from the gas-side can be 
predicted with a confidence of about 40% at moderate acoustic Reynolds numbers. 

© 2012 Elsevier Ltd. All rights reserved. 


1. Introduction 

Thermoacoustic technologies are concerned with developing new 
concepts of engines, coolers and heat pumps which operate on the 
basis of a range of thermoacoustic effects. These are broadly under¬ 
stood as energy transfer between a compressible fluid and solid 
boundary in the presence of an acoustic wave. The correct phasing 
between the acoustically-induced fluid displacement and its 
compression or expansion, coupled with heat transfer processes 
between the solid and fluid, will lead to the implementation of 
Stirling-like thermodynamic cycles which are of practical engineering 
importance. A number of practical thermoacoustic devices, both 
standing- and travelling-wave, have already been demonstrated [1,2]. 

There are significant advantages of utilising thermoacoustic 
technologies, in certain areas of application. Firstly, as there are no 
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moving components, the thermoacoustic devices can be around 
one order of magnitude cheaper than many conventional cycles 
(e.g. typical Stirling machine) and can offer longevity and low 
maintenance costs. In addition, as the thermodynamic cycle uses 
environmentally benign noble/inert gases as the working medium, 
there are significant environmental benefits of the technology. 
Finally, thermoacoustic systems scale very well and are ideal for 
utilising the low grade heat input (waste heat from industrial 
processes and solar energy and other forms of energy harvesting) 
which offers new opportunities for energy conservation. 

A typical thermoacoustic device consists of (i) an acoustic 
network, (ii) an electro-acoustic transducer, (iii) a porous solid 
medium (namely a regenerator in travelling-wave systems [3] or 
a stack in standing-wave systems [4]) and (iv) at least a pair of heat 
exchangers (HXs) [5]. The stack/regenerator is the component 
where the desired heat/sound energy conversion takes place. “Hot” 
and “cold” heat exchangers (CHX and HHX, respectively), placed in 
close proximity of both ends of the stack/regenerator, absorb (or 
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Nomenclature 

X 

axial coordinate/direction (m) 




HXs coordinate (m) 

a 

speed of sound (m s -1 ) 

y 

transverse coordinate/direction (m) 

A 

surface area (m 2 ) 

yo 

half fin spacing (m) 

Cp 

isobaric specific heat (J kg _1 K _1 ) 



d 

dimension of the square side of the cross section of the 

Greek symbols 


resonator 

a 

thermal expansion coefficient (K _1 ) 

e 

energy flux density (Wm‘ 2 ) 

T 

ratio of isobaric to isochoric specific heat 

e 

time-averaged energy flux density (Witt 2 ) 

5 

penetration depth, m 

f 

spatially averaged h function, resonance frequency 

A z 

HX length along the z direction 


(Hz) 

V 

shear viscosity coefficient (N m 2 s -1 ) 

h 

thermoviscous function, specific enthalpy (J m 3 ), heat 

6 

dimensionless temperature, (T - T W )/(T H - T c ) 


transfer coefficient (Wm 2 I< -1 ) 

K 

gas thermal diffusivity (m 2 s -1 ) 

j 

imaginary unit 

X 

wavelength of the sound wave (m) 

k 

wave number (m -1 ), thermoviscous function (m) 

V 

kinematic viscosity (m 2 s _1 ) 

I< 

gas thermal conductivity (Wm _1 K _1 ) 

? 

generic acoustic variable 

l 

half fin thickness (m) 

fc 

displacement amplitude at the joint and centre of the 

L 

heat exchanger length (m) 


channel (m) 

Nu 

Nusselt number 

p 

density of the gas (kg m“ 3 ) 

Pi 

acoustic pressure amplitude (Pa) 

a 

component of the viscosity stress tensor (N m -2 s -2 ) 

Pc 

pressure amplitude at the joint and centre of the 

T 

component of the viscosity stress tensor (N m -2 s -2 ) 


channel (Pa) 

0 

phase of the acoustic cycle (deg, rad) 

Pa 

pressure amplitude at a pressure anti-node (Pa) 

CO 

angular frequency (rad s -1 ) 

Pr 

Prandtl number 

Q 

blockage ratio 

Q 

time-averaged heat flux density (W m~ 2 ) 



Re 

Reynolds number 

Subscripts 

t 

time (s) 

0 

mean, time-averaged 

T 

instantaneous fluid temperature (K) 

1 

first order acoustic variable 

Ti 

amplitude of oscillating temperature (K) 

A 

acoustic amplitude at a pressure anti-node 

T w 

wall (fin) temperature (K) 

HX 

referring to the HX 

T h 

temperature of the hot fin (K) 

P 

isobaric 

Tc 

temperature of the cold fin (K) 

Res 

resonator 

T01-T12 

thermocouples 

X 

longitudinal, x-component 

u 

specific internal energy (J m“ 3 ) 

Y 

transverse, y-component 

Vxl 

amplitude of longitudinal velocity (m s _1 ) 

Z 

along the z direction 

Vyl 

amplitude of transverse velocity (m s -1 ) 

K 

thermal 

v c 

velocity amplitude at the joint and centre of the 

V 

viscous 


channel (m s -1 ) 




supply) energy from (or to) its ends thus enabling heat communi¬ 
cation with external heat sources and sinks. 

It is worth noting at this point that compared to the well-known 
types of HXs for steady flows HXs in the thermoacoustic context 
provide a range of new challenges. This is due to the fact that the 
flow is oscillatory. While the existing knowledge for steady flow 
arrangements is of little practical value [6,7], the relevant investi¬ 
gations in the oscillatory flow are very scarce [8-10]. Therefore, this 
work is driven by the need for a better understanding and quan¬ 
titative description of heat transfer under oscillatory flow condi¬ 
tions in thermoacoustic HXs with the aim of improving the overall 
system performance. The research presented in this paper had two 
types of activities carried out in tandem. Firstly, experimental 
facilities had to be developed for careful investigations of the 
thermal-fluid processes within the individual channels of parallel- 
plate HXs. Secondly, CFD modelling capabilities had to be devel¬ 
oped in order to capture the physics of thermal-fluid processes and 
to validate the models against the experimental data. In order to 
simplify as much as possible the experimental implementation and 
the numerical modelling the physical arrangement of the HXs was 
simplified compared to real thermoacoustic devices. In particular, 
only a pair of “mock-up” HXs is studied, with no stack/regenerator 
installed between them. These parallel-plate type HXs are arranged 
side by side with no gap between them. 


The experimental activities involved two-dimensional temper¬ 
ature and velocity field measurements using Planar Laser Induced 
Fluorescence (PLIF) and Particle Image Velocimetry (PIV) tech¬ 
niques, respectively, to obtain spatially and temporally resolved 
temperature and velocity fields within the thermoacoustic HX 
samples. On the basis of recorded temperature fields, the experi¬ 
mental data were processed to obtain the local and global, phase- 
dependent heat transfer rates and Nusselt numbers, and their 
dependence on the Reynolds number of the oscillating flow. The 
CFD modelling entailed the development of a two-dimensional low 
Mach number computational model to analyse the time-averaged 
temperature field and heat transfer rates in a representative 
domain of the HXs. These were generated by integrating the ther¬ 
moacoustic equations of the standard linear theory into an energy 
balance-based numerical calculus scheme. The model accounts for 
hydrodynamic heat transport along the transverse direction normal 
to the fin surfaces, temperature dependent thermophysical gas/ 
solid parameters and flow losses. The comparison between simu¬ 
lation results and experiments is based on the phase-averaged 
temperature fields and on the integrated data of heat flux and 
Nusselt number. This work is the first attempt of its kind to validate 
the heat transfer prediction codes in the context of a thermoa¬ 
coustic device by the detailed in-situ measurements of temperature 
fields obtained for an individual pore/channel of the heat exchanger 
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assembly. Therefore, it is hoped that the results of this work 
constitute an important scientific and technological foundation for 
future improvements of the performance of thermoacoustic 
engines where this is dependent on heat communication. 

2. Experimental setup 

The general schematic of the experimental apparatus and the PLIF 
instrumentation is shown in Fig. 1. The experimental apparatus 
consists of a 7.4 m long resonator of a square internal cross section. It 
is joined with a cubical loudspeaker box with a side length of 0.6 m 
through a 0.3 m long transition part, which increases the cross section 
from 134 mm x 134 mm on its resonator end to 500 mm x 500 mm 
on the box end. The resonator is filled with nitrogen at atmospheric 
pressure and room temperature. A subwoofer type loudspeaker is 
used to generate a standing wave of one quarter wavelength mode in 
the resonator, at a frequency (f) of 13.1 Hz. 

A sinusoidal signal from a function generator (TTi TG1010A) is 
sent to a power amplifier, whose output is then connected to the 
loudspeaker. The acoustic oscillation amplitude is varied by 
controlling the amplitude of the sinusoidal signal from the function 
generator. The acoustic pressure is monitored by a microphone 
(B&K Model 4136) flush-mounted on the end plate of the resonator. 
The pressure signal from the microphone is also used as a phase 
reference for Planar Laser Induced Fluorescence (PLIF) imaging. 

A pair of “mock-up” HXs was placed in the resonator, about 
4.6 m (0.17 of the wavelength) from the closed end of the resonator. 
The configuration of the HX pair is schematically shown in Fig. 2a, 
with a photograph of the actual implementation shown in Fig. 2b. 
Each HX consists of five active fins and eight dummy aluminium 
fins to fill up the whole cross section of the resonator. Each hot HX 
fin is made out of a brass plate with a cable heater embedded in 
a prepared channel (Fig. 2c). High-temperature epoxy was applied 
to fill up the channel and form a flat surface. Each cold HX fin is also 
made out of a brass plate, with a wide meandering channel 
prepared for water circulation (Fig. 2c). Water inlet and outlet are 
connected to the channel ends via steel tubing. The channel was 
covered with a thin brass sheet, which was soldered on top; the 
whole assembly was further machined to obtain a flat surface. All 
fins are 3.2 mm thick (2/) and 35 mm long (L) in the direction of the 
oscillating flow, with a span-wise width (Az) of 132 mm. 

Ceramic spacers (6 mm x 10 mm x 50 mm) were used in the 
front, as well as in the back (not shown) to maintain the channel 
height (2y 0 ) at 6.0 mm, and also to prevent the channel from severe 
deformation due to the thermal stresses in the fins. These ceramic 
blocks and HX fins were strengthened by supporting pillars. Only 
one channel, in the centre of the resonator, was left unobstructed 
for PLIF imaging. 

The surface temperature of the HX fins that form the channel 
under investigation were monitored by using twelve K-type ther¬ 
mocouples (Fig. 2a). They are mounted in the centre of the fins in the 



Fig. 1 . Schematic drawing of the experimental apparatus with a PLIF system. 


span-wise direction. Four thermocouples (i.e. “T01”, “T04”, “T09” 
and “T12”) are 3 mm away from the joint between hot and cold fins. 
Another four are in the middle of the cold (“T02” and “T05”) and hot 
(“T08” and “Til”) fins. The last four (“T03”, “T06”, “T07” and “T10”) 
are placed 3 mm away from the far edge of the cold and hot fins. The 
heating on the hot fins was controlled by a PID-mode temperature 
controller, which maintained a preset temperature value (e.g. 
200 °C) at the monitored point (“T07”) within the range of ±1 °C. 
Due to the thermal contact on the joints between hot and cold fins in 
the present arrangement, unwanted heat conduction was intro¬ 
duced, which makes the surface temperatures on the hot and cold 
fins decrease from left to right. The maximum temperature differ¬ 
ence between “T07” and “T09” or between “T10” and “T12” can be 
up to 20 °C. It is known from tests that the surface temperature on 
the fins along the span-wise direction was nearly constant. 

Before the temperature field measurement, the flow velocity 
fields were obtained with Particle Image Velocimetry (PIV) tech¬ 
niques. The velocity amplitude (v c ), at the joint and centre of the 
channel, and the corresponding displacement amplitude of gas 
parcel (Jc = vc/w, where co is the angular frequency) were varied 
from 1.3 m/s to 3.84 m/s and from 16 mm to 46 mm respectively. 
The corresponding acoustic Reynolds numbers based on the 
thickness (21) of the fin (Re = Vc2 l/v) -v being the kinematic vis¬ 
cosity-varied from 250 to 702. 

Temperature field measurements were then carried out at the 
same oscillation amplitude, by employing a PLIF system (LaVision) 
as shown in Fig. 1. Acetone was used as the fluorescent tracer. UV 
light sheet at 266 nm is generated by a Nd:YAG laser. It enters the 
resonator perpendicularly to its axis through a quartz window of 
3 mm thickness, is reflected by a UV mirror and becomes parallel to 
the resonator axis and normal to the surface of the heat exchanger 
plates. The laser sheet plane was 33 mm away from the front end of 
the heat exchangers, rather than at the centre of the HXs (Fig. 1 ), to 
avoid the effect of the surface mounted thermocouples on the flow 
and temperature fields. Two UV mirrors were used, to enable the 
investigation of the temperature fields inside and outside of the 
channel. The laser could be moved along the resonator side wall by 
using an automated traverse system to use either of the two UV 
mirrors as necessary. The UV mirrors are placed 300 mm away from 
the heat exchanger not to disturb the flow and temperature fields 
around the heat exchangers. Quartz windows are mounted on the 
resonator at the locations of the heat exchangers and the mirrors 
for their higher transmission to UV and fluorescence emissions. 

Acetone of optimal concentration was seeded into the flow by 
a purpose-designed seeding mechanism. PLIF images were 
captured by an intensified CCD camera. Images with effective pixels 
of 1024 x 1024 were obtained, giving a spatial resolution of 64 pm 
per pixel. Images were post-processed using LaVision DaVis 7.2 
software package. Three views were recorded with one next to 
another in the longitudinal direction to cover the area stretching 
beyond the length of the channel in front of the hot/cold HX fins. 
Prior to the acquisition of experimental measurement of temper¬ 
ature fields, additional calibration procedure was undertaken to 
establish the correlation between the fluorescence intensity and 
the fluid temperature. Further information about the temperature 
measurement with PLIF techniques can be found in [11]. 

Measurements of temperature and velocity fields were carried 
out for 20 phases within an acoustic cycle as illustrated in Fig. 3, 
which shows the measured reference pressure pc, the measured 
velocity vc and the calculated displacement fc- The directions of the 
velocity and displacement are defined within the coordinate 
system given in Fig. 1. For each phase, 100 frames were taken to 
calculate the phase-averaged temperature field. The increase of the 
number of frames beyond 100 provides negligible accuracy gains in 
terms of the resulting temperature fields. 
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3.2 mm 
6.0 mm 


c 



Fig. 2. Schematic of the configuration of the heat exchangers (a); photograph of the physical implementation (b); photograph of machined hot (left) and cold (right) heat exchanger 
fins (c). 


3. Computational model 

The numerical model used in this study is the same as that 
applied in [12], so a more synthetic description is given here. For 
additional information the reader is directed to [12]. The periodicity 
of the HXs structure along the transverse y direction allows 
calculations to be performed in a single channel of the HX 
assembly. The computational domain is further reduced by 
symmetry from half a gas duct to the solid surfaces of the coupled 
hot and cold fins which are treated as two isothermal surfaces 
respectively at temperature T H and T c (see Fig. 4). 

The implementation of the numerical model relies on a number 
of assumptions: The problem is considered as two-dimensional. 



Phase reference, deg 


Fig. 3. Reference pressure, displacement and velocity variation in an acoustic cycle and 
the investigated 20 phases. The gas displacement amplitude is 16 mm. 


The processes taking place are within a low Mach number 
regime, so that any acoustic variable J can be expressed in the 
complex notation by the conventional first-order expansion 
?(x,y, t) = ? 0 (x,y) +fte{?i(x,y)e Jwt }, where t is the time, j the 
imaginary unit, w the angular frequency, Jo the (real) mean value, Ji 
the (complex) amplitude of oscillation, and where Re {} denotes the 
real part. The fluid is Newtonian and obeys the ideal gas state 
equation. Temperature dependent thermophysical gas parameters 
are used in the model. FIXs are significantly shorter than the 
acoustic wavelength and acoustically non-intrusive so that the 
acoustic field can be approximated in the HX neighbourhood by 
a ID lossless standing wave: 

Pi = P A C0S(l<x s ) = Po V X] = j^sm(kxs) = jv 0 (1) 

where p\ is the amplitude of the dynamic pressure, Pa is the 
amplitude of the dynamic pressure at a pressure anti-node, v x \ is 
the amplitude of the longitudinal particle velocity, k is the wave 



computation domain symmetry line 


Fig. 4. Magnified view of the computational domain (light grey area). 
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number (/< = 2 tu/A), a is the sound velocity and x s the distance of the 
HXs from the closed end of the resonator. 

Amplitudes of the oscillating temperature T\, longitudinal 
velocity v x \, transverse gradient of the longitudinal velocity 0 v x i/ 0 y, 
transverse velocity v y \ and transverse gradient of the transverse 
velocity 0 v y i/ 0 j/ inside a gas pore are described by the standard 
equations of the classical thermoacoustic theory. According to the 
assumption that the specific heats of the plate/fin materials are 
notably greater than the isobaric specific heat of the gas c p (i.e. 
temperature oscillations in the solid are negligible), the following 
equations can be formulated [ 1 , 6 ]: 


1 


Po c p 
- Pr(l 


(1 -h K )p, - 


1 


dpi ar 0 


p 0 w 2 (l -Pr) dx ax 


[(i - M 


ft»)] 



v x \ 


j dp! 

cop 0 dx 


(1 




0/1 

at U 


pv z + pu) = -Ve 


(ii) 


where u is the specific internal energy and where e, the energy flux 
density, is defined as 


e - pv[ iy 2 + h ) v ct KVT 


( 12 ) 


where h is the specific enthalpy and a the viscous stress tensor. 
Time averaging Eq. (10) over an acoustic cycle leads to 


V-e = 0 (13) 

where over-bar means time-averaged. Eq. (13) is applied to each 
cell of the computational domain to impose the local energy 
balance. To accomplish this, a finite difference technique is 
employed where the quantitative results of standard linear theory 
are substituted for the x and y the components of the time- 
averaged energy flux density [13]. 




cop 0 d 


2 dgi 

2 dx 

V 


k 


V 


(4) 


. (o f [1 + (y — 1)/ K ](y — k v ) — [y + (y — l)fc K ](l —/v) 

V y1 =J—3 

Po a 


Pi 


(1 ~fv) 

■ /I f/v(y — k<) —/ic(y — kv) + (^k - kv) 1 9T 0 dpi 

+j p 0 w(l-Pr)\ (l-/v) J 9x dx 

(5) 


9fyi (o ( [1 + (y — 1 Ml — /iy) — [1 + (y — l)k K ](l — fv) 


9y Pod 1 


+] 


P 


p 0 w(l -Pr) _ 


(1 -fv) 

7v(l —^k) — /k(1 — hv) + (/Ik — kv) 


(l-/v) 


dp-i 
9x dx 


( 6 ) 


where 


cosh[(l +j)y/d K ] , 
cosh[(l +j)y 0 /<5 K ] 


cosh[(l +j)y/5 v ] 
cosh[(l +j)y 0 /<5 v ] 


(7) 


_ 8 k sinh[(l +j)y/8 K ] _ 5„ sinh[(l +j)y/<5 v ] 

/vk; — ——;—_i_ r / -i—;—777 77—T /vy — 


1 +jcosh[(l +j)y 0 /<5 K ] 


1 +jcosh[(l +j)y 0 /8 v ] 


(8) 


/ 2tz/m 2tt/w 

J pvx (^V 2 + d t — J (vx&xx + Vy^xy) 

\ 0 0 


e x 


CO 

2tu 


2tz/oj \ 

xd t-K I ^d t 


Ox 


o 


/ 


/ 2tt/w 


- CO 

6y = 2tt 


1 


2tt/w 


pVy | —v + h | dt 


V 0 


/ <r* T *y + Wyy) 


0 


2tt/w \ 

xd t-K I ^d t 


o 


9y 


/ 


where 


2 /- 0l7 x 0l7y 

ffj “ = 3n 2 ^—^ 


2 dv 


-nV 


y 


3' 0iy 


Ow 0 f v \ dv x 
Txy — Tyx — V [ 7777 T- 7777 ) ~ V 


dy dx 


9y 


(14) 


(15) 


(16) 


(17) 


_ tanh[(l +j)y 0 /k] f _ tanh[(l +j)y 0 /<5 v ] (q) 

jK [(i +j)y 0 /4] Jv [(i +j)y 0 /<5v] 1 J 

and where Pr is the Prandtl number, (3 the thermal expansion 
coefficient, y the ratio of isobaric to isochoric specific heat, po the 
mean density, To the time-averaged temperature, d v = y/2p/p^co 
the viscous penetration depth (77 being the shear viscosity coeffi¬ 
cient) and <5 k = y/2I</p 0 c P co the thermal penetration depth (I< 
being the gas thermal conductivity coefficient). 

The longitudinal pressure gradient dpi/dx along the HX 
assembly is calculated by imposing continuity of volume flow rate 
at the entrance of the HXs 


°yy 



017 x \ 4 0 Vy 

dx) ~3 71 dy 


(18) 


the x velocity gradients having being neglected because of the order 
of <5 V /A been smaller than the y velocity gradients. Making use of the 
complex notation and retaining only terms up to second order Eqs. 
(14) and (15) can be written as: 



(19) 


dpi Ppu 

dx (1 -f v )Q V ° 


( 10 ) 


where the blockage ratio Q = AnxMres = 1/(1 + //yo) describes the 
porosity of the HXs (A res being the cross sectional area of the 
resonator and Ah x is the cross sectional area of the HXs open to gas 
flow). The gas temperature is governed by the general equation of 
energy conservation [13]: 



= Ipo CpRe{T, v y ,} - IpRefei> xl j 



dVy, 

9y 



- 1< 


97b 

9y 


( 20 ) 


where energy flux densities have been converted to heat flux 
densities ( e x = q x ey = q y ) because in the short stack approximation 
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the work flux contribution to the total energy flux is small [1 ]. The 
explicit expressions of these equations are reported in [12]. 

The calculation of the steady-state two-dimensional time- 
averaged temperature distribution is performed using a finite 
difference methodology where temperature spatial gradients are 
discretized using first order nodal temperature differences. To this 
end, the computational domain is subdivided using a rectangular 
grid. In the x direction the computation mesh size is typically 
0.005 L while in the y direction the computation mesh size is 
typically 0.02 yo. The imposed boundary conditions took into 
account: the symmetry on the nodal line aty = 0; the continuity of 
temperature and transverse heat fluxes at the gas-solid interfaces 
(y = y 0 ); and the vanishing energy flux at the HX pore ends (x = 0; 
x = 21) facing the duct (gas oscillations outside the HXs are adia¬ 
batic and thus the thermoacoustic effect vanishes). 

Applying Eq. (13) in each cell of the computational grid, a system 
of quadratic algebraic equations with respect of the unknown 
variable T 0 is generated. The system is solved by a code developed 
by the authors in FORTRAN-90 language which executes the 
recursive Newton-Raphson method [14]. At each iteration the 
system of linear algebraic equations providing the temperature 
corrections for the next step is solved using a LU decomposition 
with partial pivoting and row interchanges matrix factorization 
routine. The latter is taken from the LAPACK library routines 
available online at [15], where details of the accuracy, computa¬ 
tional cost, etc. can be found. Once the time-averaged temperature 
distribution is known, it can be substituted in Eqs. (19) and (20) to 
determine the energy flux distributions along the x and y directions 
in the gas. In the next section the code is validated by comparing 
numerical predictions against experimental measurements. 

4. Results and discussion 

The temperature field measurements were performed in 20 
phases within an acoustic cycle. Here, the most characteristic 
temperature distributions at six phases are presented to discuss the 
heat transfer process in an acoustic cycle. Fig. 5 illustrates the 
temperature distribution in these six phases, $i(0°), $6(90°), 
$s(126°), $12(198°), $16(270°) and $i S (306°), for a Reynolds 
number equal to 250 (with the corresponding fluid displacement 
amplitude of 16 mm). 

At phase $i, the fluid displacement reaches the maximum 
negative value, which indicates that the gas parcel moves the 
farthest to the left. The cold gas in the channel of the cold fin 
(referred as the cold channel thereafter) moves far into the channel 
of the hot fin (referred as the hot channel thereafter) and is heated 
by the plates. Some part of hot gas in the hot channel flows out of it 
and enters the open area behind the hot fins. It can be clearly seen 
that the gas in the top region behind the hot fins is heated up, due to 
the cumulative effect of natural convection over time. After this 
phase, the gas parcel starts to accelerate to move to the right. At 
phase $6, the cold gas heated by the hot fins partly flows out of the 
hot channel and enters into the cold channel. Hot gas entering from 
the open area behind the hot fins is dominant at the left end of hot 
channel. The gas parcel continues to move to the right. At phase $g, 
the hot gas has dominated the whole hot channel and partly 
penetrates into the cold channel. The cold gas in the cold channel 
partly moves into the open area behind the cold fins and mixes with 
the colder gas in this region. It can be found that the natural 
convection effect also exists. At phase $ 12 , the gas parcel just starts 
to move to the left. The cold gas in the cold channel gradually moves 
into the hot channel and, at the same time, the cold gas in the open 
area behind the cold fins partly flows into the cold channel. At phase 
$i 6 , the cold gas penetrates far into the hot channel and the hot gas 
in the hot channel also migrates into the open area behind the hot 
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Fig. 5. Temperature distributions for six selected phases and for Re = 250. Red rect¬ 
angles mark the heated plates and blue rectangles mark the cooled plates in the heat 
exchanger assembly. The black dashed lines mark the position where the three views 
(left, central and right) join to one another. (For interpretation of the references to 
colour in this figure legend, the reader is referred to the web version of this article.) 


fins. At phase $ig the hot channel is occupied by the cold gas and the 
hot gas in the hot channel moves further to the left. After this phase, 
the gas parcel keeps moving to the left. Then another acoustic cycle 
repeats itself. Analogous processes can be identified at the higher 
Reynolds number where the relatively large velocity develops 
a very thin thermal boundary layer. 

For comparison, in Fig. 6 the time-averaged temperature 
distribution generated by the computational model is shown for 
the same value of the Reynolds number. At the joint position and in 
the mid cannel regions (far from the fin surface) the temperature is 
nearly uniform over a cross section (isotherms are vertical). In the 
direction away from the joint towards both the cold and hot 
channel the isotherms become even more curved and practically 
horizontal at the gas-solid interface. Isotherms, however, become 
more densely packed when approaching the joint point and this 
trend is compatible with a growing transverse temperature 
gradient dT 0 /dy near this region. The joint involves a great 
temperature gradient in both the longitudinal and the transverse 
direction and most of the heat transfer between the gas and the 
plate should be localized in this region. This is confirmed by the 
longitudinal distribution of the transverse heat flux density at the 
gas-solid interface calculated by Eq. (19). Results are shown in 
Fig. 7 at selected Reynolds numbers. 

Since the time-averaged transverse energy flux at the plate 
surface is positive when entering the plate and negative when 
leaving the plate, Fig. 7 implies that over the period of an acoustic 
cycle, energy flows out of the hot fin on the left, then hydrody- 
namically along the thermal boundary layer in the gas and into the 
cold fin on the right. 

Information on the time evolution of the heat flux exchanged 
between the oscillating gas and the solid fins can be obtained from 
the cross sectional temperature profiles. In Fig. 8 the transverse 
temperature profiles measured 1 mm away from the centre of the 
HX assembly both in the cold and hot fin are illustrated for the 
previously selected six phases and for Re = 250. The six phases 
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Fig. 6. Time averaged temperature distribution in half a gas channel generated by the 
computational model for Re = 250. 
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denoted with red symbols correspond to the position in the hot 
channel 1 mm away from the “joint”, while the other six phases 
denoted with blue symbols correspond to the position in the cold 
channel 1 mm away from the “joint”. The distance from the fin 
surface is normalized by the channel width, 2yo = 6 mm while the 
gas temperatures are normalized by 


0 


T-T w 

Th-T c 


( 21 ) 


T, T w , 7 h and Tc being the instantaneous gas temperature, wall (fin) 
temperature, reference hot temperature (200 °C) and reference 
cold temperature (30 °C), respectively. The temperature profiles 
clearly show the large temperature variations in a very thin layer of 
gas near the wall. For all six phases in the hot channel the heat 
transfer direction is from the plate to the gas, while in the cold 
channel the heat is always transferred from the gas to the plate. It is 
also found that in the cold channel at phase $i, the temperature of 
the gas beyond the thin layer near the wall is lower than the wall 
temperature (a negative 0). It is understood that at this phase the 
gas parcel is in the left-most position and the cold gas from the cold 
channel moves deeply into the hot channel. The colder gas from the 
right takes over this position. This part of gas affects the processes 
for a while so that at phase $6 the temperature in the centre of the 
cold channel is still lower than the wall temperature. 

The temperature profiles enable calculation of local, phase- 
dependent, transverse heat flux densities on the HX surface from 
first principles: 



Fig. 7. The time-averaged transverse heat flux density at the hot and cold fin surface 
(y = y 0 ) as a function of the axial coordinate for selected Reynolds numbers. 


Fig. 8. Cross sectional temperature profiles measured 1 mm away from center of the 
HX assembly both in the cold and in the hot fin for the previous selected six phases and 
for Re = 250. 


q y (x, <S>) 


d T(x,y,<Z>) 

dy 


y=yo 


( 22 ) 


Integrating Eq. (22) in phase over an acoustic cycle and/or in 
space over the fin length it is possible to obtain time-averaged and/ 
or space-averaged heat flux densities. In particular, phase-and-area 
averaged heat transfer coefficients are calculated as 




Qy(x i T>)/ (T w 


T) 


d<2>. 


(23) 


T is the reference fluid temperature evaluated in the centre of 
the channel, at the same stream-wise location as the joint between 
hot and cold channels. The values of h can then be converted into 
non-dimensional Nusselt numbers based on the fin thickness 21: 



(24) 


Note that when manipulating data generated by the computa¬ 
tional model the phase integration is not required since the model 
directly provides time-averaged quantities. 

The comparison between the experimental and computational 
data is presented in Fig. 9 in terms of the area-and-phase averaged 
transverse heat flux density at the fin surface. The graph clearly 
shows that the two fins exhibit a rather different behaviour. The heat 
flux exchanged at the hot HX (open circles) appears to be an 
increasing function of the Reynolds number. This is due to the high 
velocity and displacement amplitudes at high Reynolds numbers 
which increase the temperature gradient in the hot and cold 
channel. The heat flux exchanged by the cold HX (full circles) rea¬ 
ches a maximum for Re = 535. This behaviour can be interpreted as 
an effect of the variation of the effective length of heat transfer with 
the gas displacement amplitude. At low Reynolds number, in fact, 
the gas displacement amplitude is less than the length of hot and 
cold HX, so the gas portion heated by the hot exchanger can only 
partly transfer heat to the cold HX. In particular, the gas heated at the 
left-most end of the hot HX cannot exchange heat with the cold HX. 

When the Reynolds number is 535, the corresponding gas 
displacement amplitude is 36 mm, which is close to the length of 
hot and cold HX (35 mm). In this case the gas heated by the hot HX 
can fully transfer heat to the cold HX. When the gas displacement 
exceeds the length of hot and cold heat HXs, the colder gas from 
outside of the hot HX moves into the cold HX and cannot effectively 
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Fig. 9. Area-and-phase averaged transverse heat flux at the cold (full circles) and hot 
(open circles) heat exchanger. Solid lines are predictions of the simulation model. 



Fig. 10. Area-and-phase averaged transverse heat flux at the cold (full circles) and hot 
(open circles) HX. Solid lines are predictions of the simulation model modified to 
account for heat leakage to the resonator duct. 


exchange heat with the cold HX. In the considered Re range the 
experimental heat fluxes are matched by the computational ones 
with an error of about 40%. The model is able to take into account 
the “saturation” of the heat transfer at increasing Re numbers (the 
growth of the curves diminishes with Re) but not the maximum 
observed in q y at the cold HX. Furthermore, in contrast with 
experiment, the heat exchanged by the cold HX appears to be 
higher than the heat exchanged by the hot HX. This result is not 
surprising since in the proposed model the channel ends are 
considered as thermally insulated, based on the hypothesis that 
oscillations in fluid outside the HXs are adiabatic, and so no energy 
can escape from the channel. Thus, the cold HX is loaded by the 
heat delivered by the hot HX plus eventually the heat generated by 
viscous dissipation at the solid walls. 

However, Fig. 5 clearly shows that a large fraction of the heat 
transmitted to the gas by the hot HX is spread out in the resonator 
hot duct by free convection. This phenomenon is not negligible at 
the hot duct side since its strength is proportional to d 3 , 
(d = 134 mm being the side of the resonator’s square cross section) 
and to the difference between the mean temperature of the gas 
leaving the hot HX pores and the mean temperature of the hot duct 
(approximately the difference between 150 and 200 °C). At the cold 
duct side the free convection process is much less marked due to 
the lower temperature difference between the gas and cold duct 
(nearly zero). Inside the HX channels the phenomenon is clearly 
negligible since 2yo << d, as confirmed by inspection of Fig. 5 that 
reveals that the gas-fin heat exchange is symmetric along the 
transverse y direction. This provides an “a posteriori" justification of 
the reduced computational domain selected on the basis of the 
(assumed) symmetry of the problem along the y coordinate. 

The strong free convection processes taking place at the hot duct 
side facing the HXs assembly entail, however, that the hypothesis of 
considering the gas oscillations just outside the hot HX fins as 
adiabatic (applied in the implementation of the model) is probably 
no longer valid. These circumstances are responsible for the higher 
measured heat transfer rates at the hot HX compared to those 
computed by the model, and for the breakdown of the y symmetry 
at the inlet and outlet of the HX pores. 

Therefore, as a first step to take into account these very complex 
free convection processes in the framework of the simplified 
computational model used, the numerical code has been modified 
by introducing two fictitious thermal reservoirs on both the side of 
the cold duct and the side of the hot duct. These thermal reservoirs 


are modelled to thermally interact with the oscillating gas inside 
the HX channel through two heat transfer coefficients that, in 
addition to the thermal reservoir temperatures, are considered as 
fitting parameters. Results are shown in Figs. 10 and 11 in terms of 
the transverse heat fluxes and the Nusselt numbers. The corre¬ 
sponding heat leakages to the hot and cold side of the duct 
computed by the simulation model are reported in Fig. 12. 

The experimental data are now described by the simulation 
curves with an error of 3% and 35% at the hot and cold HX, 
respectively. The measured Nusselt numbers exhibit a dependence 
on Re similar to that of q y (Fig. 11). A maximum in Nu evaluated at 
the cold HX for Re = 535 is observed to which a saturation in Nu 
evaluated at the hot HX corresponds. The experimental Nu 
numbers are fitted by the numerical model with an error of 9% and 
35% at the hot and cold HX, respectively. 

The previous analysis suggests that in the investigated range of 
flow conditions, Re = 535 is the optimal flow parameter for the 
thermal performance of the tested HXs. This is likely due to the 
circumstance that the gas displacement amplitude equals the 



Fig. 11. Area-and-phase averaged Nusselt numbers at the cold (full circles) and hot 
(open circles) heat exchanger. Solid lines are predictions of the simulation model 
modified to account for heat leakage to the resonator duct. 
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Fig. 12. Numerically computed heat leakages to the hot side and cold side of the 
resonator duct. 

length of the hot and cold exchanger for this Reynolds number. 
Heat transfer coefficients from the gas-side can be predicted with 
a confidence of about 40% at moderate acoustic Reynolds numbers. 
The matching of the experimental data improves if heat transport 
from the couple of HXs to the surrounding environment (hot and 
cold ducts) is considered. These heat losses are most likely caused 
by the entrance/exit effects localized at the resonator-HX cross 
section interfaces which are responsible of complex non-linear 
temperature and flow patterns (strong vortex flows and turbu¬ 
lence generation effects [16,17]). Although the simplified model 
proposed to account for this effect doesn’t constitute at this stage 
a predictive algorithm for quantitative estimations, it provides 
some evidence that this kind of phenomena could have a non- 
negligible influence on the value of the heat transfer coefficients 
at the gas side of the HXs (especially at the hot HX). Further 
experimental investigations and theoretical modelling are required 
for more detailed insight into the physics of thermo-fluid processes. 

5. Conclusions 

The work described in this paper utilises the novel measure¬ 
ment techniques developed in order to obtain instantaneous 
temperature and velocity fields in the interstices of the internal 
structures of thermoacoustic devices. This is a challenging task, 
which has not been attempted before, and which enables a unique 
insight into the physics of thermofluids processes within the 
oscillatory compressible flows that govern the associated energy 
transfer. Acetone-based Planar Laser Induced Fluorescence (PLIF) 
and Particle Image Velocimetry (PIV) techniques have been applied 
to obtain spatially and temporally resolved temperature and 
velocity fields within parallel-plate thermoacoustic heat 
exchangers. The local temperature distribution profiles within an 
acoustic cycle and space-averaged phase-dependent heat fluxes 
reveal the heat transfer variations phase by phase. The dependence 
of the space-cycle averaged heat fluxes and the Nusselt numbers on 
the Reynolds number shows that the optimal performance of heat 
exchangers can be achieved when the gas displacement amplitude 
is close to the length of hot and cold heat exchanger. 

The experimental data obtained from the measurement of the 
temperature field provide an important reference for validation of 
a CFD code implemented on the basis of the classical thermoacoustic 


linear theory. Comparisons between the detailed spatially- and 
temporally-resolved measurements provide a unique opportunity 
to test the predictive capabilities of the CFD code developed in the 
course of this research. It is shown that the code is able to predict 
heat transfer coefficients from the gas-side with a confidence of 
about 41% at moderate acoustic Reynolds numbers. However, better 
estimates can be achieved when the entrance/exit effects localized 
at the resonator-HX cross section interfaces are taken into account. 
These effects are responsible for considerable heat losses from the 
couple of heat exchangers into the surrounding environment (the 
resonator extending beyond the heat exchangers). The results of this 
study demonstrate the need for developing the new thinking for the 
design of CFD codes to encapsulate the nature of thermal-fluid 
interactions present in the thermoacoustic devices, and in the 
long term to encourage new paradigms and concepts for the future 
design of heat exchangers in thermoacoustic systems. 
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